Fast calculation of thermodynamic and structural 
parameters of solutions using the 3DRISM model and the 

multi-grid method. 

p I Volodymyr P. Sergiievskyi 

Max Planck Institute for Mathematics in the Sciences, 
Inselstrasse 22, DE-04103, Leipzig, Germany 
e-mail: sergiiev@mis.mpg.de 

May 1, 2012 



Note: the information in this preprint is not up to date. Since the first pubhcation of the preprint (9 
Nov 2011) the algorithm was several times modified which allowed to achieve much better computational 
results. The paper which describes the new algorithm and new computational tests was recently accepted 
for the publication in the Journal of Chemical Theory and Computation as: 
Volodymyr P. Sergiievskyi, Maxim V. Fedorov 

"3DRISM Multi-grid Algorithm for Fast Solvation Free Energy Calculations" 
( DOI: 10.1021/ct200815v , http://pubs.acs.org/doi/abs/10.1021/ct200815v, ) 

Abstract 

In the paper a new method to solve the tree-dimensional reference interaction site model 
(3DRISM) integral equations is proposed. The algorithm uses the multi-grid technique which 
allows to decrease the computational expanses. 3DRISM calculations for aqueous solutions 
of four compounds (argon, water, methane, methanol) on the different grids are performed in 
order to determine a dependence of the computational error on the parameters of the grid. 
It is shown that calculations on the grid with the step O.OSA and buffer SA give the error of 
solvation free energy calculations less than 0.3 kcal/mol which is comparable to the accuracy 
of the experimental measurements. The performance of the algorithm is tested. It is shown 
that the proposed algorithm is in average more than 12 times faster than the standard Picard 
direct iteration method. 

1 Introduction 



Integral equation theory of liquids (lETL) is an effective method for the prediction of the structural 
and thermodynamic parameters of liquid and amorphous matter IHE]. lETL describes the stucture 
of the liquid by using the correlation functions. The main equation of lETL is the Ornstein- 
Zernike (OZ) equation which for the case of molecular liquids connects the correlation functions 
of six independent variables [1]. Because of the comutational complexity solving the OZ equation 
in a general case is still an open task. In practice one uses approximate models, such as the 
reference interaction site model (RISM) [3] and three-dimensional reference interaction site model 
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(3DRISM) [4J. It was shown that the RISM equation as the imphcit solvent model gives more 
physically-correct description of the system in comparison to the continuum electrostatics models 
O ini [7]. The important application of the RISM equations is the calculation of the solvation free 
energy (SFE) of compounds. Within the framework of the proposed model solvation free energy 
can be analytically calculated from the solutions of the RISM equations IH [9] . Recently there were 
proposed several parameterization methods of the results of the RISM calculations which predict 
the solvation free energy with the accuracy of 1 kcal/mol [ini ttH [ISl [131 E]- However, despite 
the good results RISM still describes not enough accurately the molecular solvent, which leads 
to the numerous corrections needed for achieving the desired accuracy of the calculations. The 
3DRISM equations describe more accurately e the stucture of the molecule and, as it was shown, 
they need essentially less corrections for accurate calculation of the thermodynamic parameters 
p31 116] . Direct Picard iteration method is the standard algorithm for solving the equations of the 
lETL [T7]. Although it is easy to implement, but it has quite a low convergence rate. There were 
proposed several alternative schemes which increase the convergence rate [H UHl [191 [20] . Another 
way to speed up the calculations is to use the multi-scale methods [211 [221 IMl [211 l25l [211 123 128] . 
It is necessary to note that although most of the multi-scale methods use several grids, not all of 
them do it effictively. However, the multi-grid technology allows to use all the advantages of the 
multi-scale approach [29] . It was recently shown that the multi-grid allows to increase performance 
of the RISM calculations up to the several dozen times [30]. In the current paper the multi-grid 
algorithm which solves the 3DRISM equations is described. The optimal discretization parameters 
for the solvation free energy calculations are found. The benchmarking of the algorithm and the 
comparison to the standard Picard iteration method is perfromed. 



2 Method 
2.1 3DRISM 

In the paper the 3DRISM for the description of the infinitely diluted solutions is used. The solvent 
molecules are given in the RISM approximation, while the solute molecule is a three-dimensional 
object. The solvent structure is described by the total and direct correlation functions of the 
solvent site a: hair), Ca(r). 3DRISM equations are written as following: 

-^solvent p 

h^{r)= J2 / c^{r')xa^' - r)dr' (1) 

where Absolvent IS the number of solvent sites, x^a{r) is the solvent susceptibility function. The 
functions X(a{r) are defined by the formula x^ai^) = uj(a{r) + ph^a{r), where h^aij') is the total 
correlation function of the pure solvent, p is the solvent number density (number of particles in 
the unit volume), u^aij) = 6^aS{r) + (1 — S^a)S{r — rgQ,)(47rr|^)~\ 6^a is the Kronecer's delta, 6{r) 
is the Dirac delta fuunction, r^a is the distance between the sites ^ and a in the solvent molecule. 
Equation ([ID is completed with the closure relation: 

Where /3 = l/Zc^T, is the Boltzmann constant, T is the temperature, Ua{r) is the interaction 
potential of the site a with the solute molecule. To solve the 3DRISM equations iteratively the 
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Figure 1: The spacing is the distance between the grid points, and the buffer is the minimal 
distance from the solute atoms to the grid boundaries 



buffer 
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equations are written in the following form [31] : 

-^solvent 



" solvent p 

7.(r)= / C[^^{r'-r)].xa^')dr' + W-Cb.{r)] 



(3) 



where 7„(r) = /i„(r)cf(r), cf(r) = c«(r) + /3t/i'(r), f/«(r) = f/f(r) + U^{y), f/f(r) is the short- 
term component of the potential, U^{y) is the long-term component of the potential, (r) = 
—(5 /jg3 U^{y — r')x^a{y')dT' In the current paper the Kovalenko-Hirata closure is used [17j. It 
is defined by the following formula: 

C[7<,(r)l = { '^sIX"' ~ " ' " ^ " " (") 

In the numerical formulation of the problem ([I])-(I2]) the functions 7a (r), ^0(1") are given 

by their values in the grid nodes. We use two values to define the grid: spacing and buffer. The 
spacing is the distance between the grid points, and the buffer is the minimal distance from the 
solute atoms to the grid boundaries (see Figure [1]) 

We use the notations -^j-^f-] for the forward and inverse discrete Fourier transforms on 

the grid Q. The discrete equation which corresponds to the equations (I3j) is the following: 



x-j-g[c[r^]] +0^-c[r^] (5) 



where = (7? , . . . , tL^)^ = {Of, • • • , ^L™„J^' = [4oW-e.x^.,_, upper index 
Q means that the functions are defined by their values in the points of the grid Q. The euqation 
([5]) can be written more compactly: 

pe = F[T^] (6) 
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where F[T^] = J^g^ X-J^g[C [T^]] + 0^ - C [r^] . The Picard iteration method is defined by 
the foUwing reccurent formula: 

r^+i = (1 - A)r^ + AF[r^] (7) 

where - is the n-th step solution approximation, A is the mixing parameter. 
2.2 The multi-grid method 

The multi-grid method is used in the paper to reduce the comutational time [291 [30]. The problem 
(|3]) is discretized on several different grids. We say that the grids which have the smaller number 
of points the "coarse" grids, and the grids which have larger number of points are the "fine" grids. 
We introduce the operators p[-], r[-] which convert the coarse grid to the finer one and the finer grid 
to the coarser one correspondingly. We introduce the operator R[-] which transforms the solution 
from the fine grid Q to the coarse grid r[Q]. 

R[T^] = r^i^i (8) 

We introduce also the operator P[-], which interpolates the coarse-grid solution to the fine grid 

g. 

p[r^[^]] = (9) 

In our work the linear interpolation operator P[-] is used. The multi-grid method solves the 
following problem: 

r^ = F^[r^] + D^ (10) 



where D = [df, . . . , ^ J is the vector of corrections. 

It is convenient to introduce the iterative operator Ag[-; ■], which is defined in the following 
way: 

Ag[Tg; D^] = (1 - A)r^ + A {Fg[T^] + D^) (11) 

The iterative multi-grid algorithm which solves the equation ffTOj) can be written in a following 
way: 

r^+i = HKDT (12) 

where is the solution approximation on the n-th step, A^g[-; ■] is the multi-grid operator, which 
performs one multi-grid iteration step of depth / on the grid Q. Calculation of the multi-grid 
operator of depth Z = is equal to the performing of n one-grid iteration steps on the grid Q If 
/ > for the given initial guess and correctons vector the multi-grid operator A^g[-; ■] is 
calculated by the following algorithm: 

Input: F^, D^, / 

Output: F^+i = M[^[F^;De] 

1. Perform Ui Picard iteration steps on the fine grid (in our work ui = 1): 

F'^ = (A,)- [F^;DT 

2. Move to the coarse grid r[Q]: 
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3. Calculate the coarse grid correction: 

— p rrrr'Gii _ PIT 

(0) 



E^i^i = R [F[r'^]] - F[rf^' 



4. Repeat recoursively /i multi-grid iteration steps of depth / — 1 on the coarse grid (in our work 
/x=l): 

5. Correct the fine-grid solution: 

F'/G _ r'G _i_ p IrriG] _ -pr[g] 

i - i + ^ [i (^) 1 (0) 

6. Perform 1/2 Picard iteration steps on the fine grid (in our work U2 = 0): 
The iterations are performed until the following stop-condition satisfies: 

llr^ r„+m|| < £tres (13) 

where m is selected in such a way that the following relation holds: 

\\-pG _ ll^nmlir^ — ll 

In the paper is used the norm based on the solvation free energy calculations: 

||r? - r^ll = lAGKHiT^) - AGKHiT^)] (15) 
where the solvation free energy is calculated using the Kovalenko-Hirata formula [32] : 

-^solvent /» -| 

AGKHir^) = pknT V / 9{-h^{r))h^{r) - -c^{r)h^{r) - c^{r)dr (16) 
In the paper the value theshold value £:tres=0.001 is used. 
2.3 Computational details 

For the argon, methane, metanol and water the OPLS-AA partial charges were used[33]. For the 
water (solute and solvent species) the parameters of the MSPC-E water model were used [M] . 
In the paper the total correlation functions of the pure water calculated in the paper |3l] were 
used. Pairwise a Lennard- Jones parameters were calculated as the arithmetic mean of the atomic 
parameters, pairwise e Lennard- Jones parameters were calculated as the geometric mean of the 
atomic parameters: 

CTu = — - — £^12 = V^i • £2 (17) 

Calculations were performed on the Dual Core AMD Opteron(tm) Processor 885 with the clocking 
2613 MHz. The forward and inverse Fourier transforms were calculated using the FFTW3 library 
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3 Results 



The 3DRISM calculations were performed for the infinitely dilutes aqueous solutions of four com- 
pounds: argon, methane, methanol, and water. To determine the optimal grid parameter the 
solvation free energy calculations for the grids with the different spacing and buffer were per- 
formed. In the Figure [2] the dependence of the calculation error on the grid spacing is presented. 
For the calculations the grids with the buffer SA and the spacing from O.OSA to 0.6A were used. 
The calculation error on each grid was calculated as the absolute value of the difference between the 
SFE calculated on the current grid and the SFE calculated on the grid with the spacing 0.025A . 
The solvation free energy values calculated on this grid were used as the standard. On the grid 
with the spacing O.OSA the maximal calculation's error is less than 0.3 kcal/mol, while for the 
grids with the larger number of points the calculation's error is more than 1 kcal/mol. Relatively 
large calculations' errors can be explained by the fact that the Cartesian grid cannot approximate 
the spherically symmetrix functions 0J^a{f) which define the structure of the solvent molecule with 
the good accuracy. Thus in a fact we use the solvent molecules with the distorted geometry which 
causes the decreasing of the accuracy of the calculations. The acceptable accuracy for the chem- 
ical applications is 0.5 kcal/mol. As it was shown, it is enough to have the grid with the spacing 
O.OSA in order to calculate the SFE with the accuracy 0.5 kcal/mol. 



Figure 2: Dependence of the error of the SFE calculations (Error of G) on the grid spacing 
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In Figure |3]the dependence of the calculation error on the grid buffer is presented. Calculations 
were performed on the grids with the constant spacing Ai?=0.05A . The error was calculated as 
the difference between the SFE calculated on the current grid and the SFE calculated in the grid 
with the spacing O.OSA and the buffer 14A which was used as the standard. For the calculations 
with the accuracy of 0.5 kcal/mol it is enough to use the buffer sA . 

To test the performance of the proposed algorithm the calculations for the aqueous solutions of 
argon, methane, methanol, and water on the grid with the spacing 0.05A and the buffer sA were 
performed. The multi-grid algorithm with the depth / = 2 was used. In Table 1 the comparison 
of the comutation time of the multi-grid and Picard iteration methods is presented. The average 
computation time for the Picard iteration method is 7 h. 42 min., while the average computation 
time for the proposed multi-grid method is 36 min. In average the multi-grid method is 12.2 times 
faster than the Picard iteration method. 
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Figure 3: Dependence of the error of the SFE calculations (Error of G) on the grid buffer. 
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Table 1: Comparison of the co mutation time for the Picard iteration method and for the multi-grid 



method. 



Compound 


Picard iteration 


multi-grid 


speedup (times) 


argon 


4 h. 45 min. 


25 min. 


11.5 


methane 


10 h. 42 min. 


41 min. 


15.7 


methanol 


12 h. 16 min. 


48 min. 


15.3 


water 


3 h. 5 min. 


29 min. 


6.4 


Average 


7 h. 42 min. 


36 min. 


12.2 
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4 Conclusions 



In the paper the new muhi-grid based method for solving the 3DRISM equations is proposed. 
To determine the optimal grid parameters the 3DRISM calculations for infinitely diluted aqueous 
solution of argon, methane, methanol and water were performed. It is shown that on the gird with 
the spacing 0.05A and the buffer SA the maximal comutational error is less than 0.3 kcal/mol, 
while for the grids with the larger spacing the computational error is more than 1 kcal/mol. The 
performance of the proposed algorithm was compared to the performance of the standard Picard 
iteration method. It is shown that the average calculation time is 36 min. for the proposed method 
while the average calculation time for the Picard iteration method is more than 7 hours. Thus, we 
show that the multi-grid algorithm is more than 12 times faster than the standard Picard method. 
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